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The tetrad-based equations for vacuum gravity published by Estabrook, Robinson, and Wahlquist 
are simplified and adapted for numerical relativity. We show that the evolution equations as par- 
tial differential equations for the Ricci rotation coefficients constitute a rather simple first-order 
symmetrizable hyperbolic system, not only for the Nester gauge condition on the acceleration and 
angular velocity of the tetrad frames considered by Estabrook et al, but also for the Lorentz gauge 
condition of van Putten and Eardley and for a fixed gauge condition. We introduce a lapse function 
and a shift vector to allow general coordinate evolution relative to the timelike congruence defined 
by the tetrad vector field. 



I. INTRODUCTION 

Gravitational wave detection is soon to become a real- 
ity. However, the scientific community has not yet been 
able to calculate gravitational waveforms from the most 
likely source, namely the violent and dynamic merger 
phase of binary black hole collisions pj. We hope to 
contribute to this effort with detailed small-scale studies. 
Our first such study of ID colliding gravitational plane 
waves indicated that two important factors which im- 
prove the accuracy and stability of the numerical calcu- 
lations are hyperbolicity of the equations and evolving 
variables which are related to physical quantities. As a 
way to generalize our results to black hole spacetimes 
and higher dimensions, we are investigating a tetrad ap- 
proach based on the formalism published by Estabrook, 
Robinson, and Wahlquist (ERW) 3J. In this paper, we 
present a modified version of the formalism, adapted for 
numerical relativity. Subsequent papers will present our 
numerical results. 

The standard approach to vacuum numerical relativ- 
ity is a 3 + 1 decomposition, of the type introduced by 
Arnowitt-Deser-Misner Q| . The 3 + 1 formulations slice 
four-dimensional spacetime into three-dimensional space- 
like hypersurfaces (see for a recent review). They 
evolve the spatial metric and extrinsic curvature of the 
hypersurfaces, which are expressed in a coordinate ba- 
sis. The evolution of the coordinates is described by a 
lapse function and shift vector, which may or may not be 
dynamic. 

A tetrad formulation uses orthonormal basis vectors, 
e a (a = 0, 1, 2, 3), which describe local Lorentz frames. 
The spacetime metric, the dot product of the basis vec- 
tors, is everywhere the Minkowski metric, g a p = e a -ep — 
rjaff . The timelike vector field of the orthonormal frames, 
eo, defines a preferred timelike congruence, to which it 
is tangent. The spatial triad vectors in a particular rota- 
tional orientation with respect to eo are e±, e-z, and e^. 
The dual basis of orthonormal one-forms is e Q such that 
(e a , ep) = Sg and e a ■ eP — rj al3 . As Estabrook and 



Wahlquist point out @, these tetrad frames are natural 
for measuring observable physical quantities. 

Most of the recent tetrad formalisms assume that the 
tetrads are tied to the physically defined flow of a fluid, 
primarily in the context of cosmology or of interior met- 
rics of rotating stars. Here, we advocate that a tetrad 
formalism may be useful even in vacuum blackhole space- 
times. The variables in a tetrad formalism are the con- 
nection coefficients (often called the Ricci rotation coeffi- 
cients), the tetrad vector components, and, typically, the 
tetrad components of the Weyl tensor or the Riemann 
tensor. The spatial coordinates are often assumed from 
the start to be comoving with the timelike congruence 
generated by the tetrad field. In dealing with black hole 
event horizons, however, it is important to allow general 
choices of coordinates. We do this by introducing a lapse 
function and a shift vector defined relative to the congru- 
ence world lines (see also van Putten and Eardley 0]). 

The ERW formulation is a quasi-FOSH (first order 
symmetric hyperbolic) system. We use the term "quasi" 
because the directional derivatives along the spatial 
tetrad legs contain partial time derivatives (also see j8]). 
The basic quasi-FOSH structure of the equations involves 
the Weyl tensor components as variables. The quasi- 
FOSH system also includes equations derived from the 
Nester gauge conditions 0. These equations evolve the 
acceleration of the congruence worldlines, a, and the an- 
gular velocity of the spatial tetrad vectors relative to 
Fermi- Walker transport, ui. The tetrad components a a 
and uj a are gauge quantities since the spacetime orienta- 
tion of the tetrad is not fixed by the spacetime geometry. 
Both a and us are orthogonal to eo, and therefore have 
only the three spatial tetrad components. In addition to 
evolution equations, the Nester gauge conditions provide 
equations which constrain the spatial dependence of the 
a a and the uj a at any given time. 

ERW's formalism can be modified to give a particularly 
simple quasi-FOSH system, as suggested by Estabrook 
and Wahlquist (EW) [10( , by eliminating the Weyl tensor 
components as separate variables, and adding the Nester 
constraint equations to the evolution equations. 
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In this paper, we discuss the derivation of the EW 
formalism in a way which is perhaps more accessible 
to those familiar with standard tensor analysis, as op- 
posed to Cartan differential form analysis. What distin- 
guishes the EW formalism is that the basic quasi-FOSH 
structure involves only the connection coefficients. Most 
of the tetrad formulations in the literature include the 
Weyl tensor or the Riemann tensor to get a simple quasi- 
hyperbolic structure. Unless one adds the constraints to 
the Einstein equations in the particular way we present 
in this paper, the quasi- hyperbolic system involving only 
the connenction coefficients is quite complicated. We do 
not claim that the EW system is simpler or more elegant 
that those involving the Weyl or Riemann tensor; how- 
ever, it has fewer variables and fewer constraint modes, 
which may be advantageous in 3D codes. 

After deriving the EW formalism, we proceed to ex- 
tend it into a form useful for vacuum numerical relativity. 
We consider two gauge conditions other than the Nester 
gauge which also give quasi-FOSH systems of equations. 
These are the Lorentz gauge used by van Putten and 
Eardley 7], and a fixed gauge, where a a and uj a are fixed 
functions of time and space. Additionally, we allow for 
a completely arbitrary relationship between the congru- 
ence and the hypersurface. Finally, we analyze the true 
hyperbolicity of the equations when expressed in terms of 
partial derivatives. The partial differential system con- 
tains evolution equations for eighteen non-gauge connec- 
tion coefficients, six gauge quantities (if a dynamic gauge 
condition is chosen), three components of a vector de- 
scribing the velocity of the congruence relative to the 
hypersurface, and nine components of the spatial tetrad 
vectors projected into a t = constant hypersurface. These 
last twelve variables, together with our lapse function and 
shift vector (which we do not evolve), completely deter- 
mine the sixteen tetrad vector components. The flux Ja- 
cobian of the partial differential system of equations has 
a complete set of eigenvectors and real eigenvalues, thus 
satisfying the requirement for hyperbolicity as per LeV- 
eque Furthermore, these equations are symmetriz- 
able hyperbolic as defined by Lindblom and Scheel 12] . 

Tetrad formulations for general relativi ty o ther than 
ERW and EW include those by Friedrich Hal . Jantzen, 
Carini, and Bini |14l. van Elst and Uggla [lJ|, Choquet- 
Bruhat and York [g, and van Putten and Eardley 0. 
Except for EW, all of these systems include the Weyl or 
Riemann tensor components as fundamental variables. 
Friedrich's paper is a definitive discussion of hyperbolic- 
ity for both tetrad and 3 + 1 representations. Janzten, 
Carini, and Bini give an extensive historical review of 
the tetrad and 3 + 1 approaches. They provide a uni- 
fied framework for the two approaches in what they call 
"gravitoelectromagnetism" , in which, however, they only 
consider co- moving spatial coordinates, van Elst and 
Uggla's formalism applies only to non-rotating congru- 
ences and orthogonal hypersurface slicings. Also, they 
do not bring the equations into a quasi-FOSH form. The 
main emphasis of Choquet-Bruhat and York's paper is 



a tetrad formulation for fluids, where the congruence is 
aligned with the fluid flow lines, and the spatial coordi- 
nates are co-moving with the congruence. The accelera- 
tion of the congruence worldlincs is given by the acceler- 
ation of the fluid, and the angular velocity variables (our 
w„'s) are fixed functions. Choquet-Bruhat and York do 
briefly mention the vacuum case, in which their system is 
quasi-FOSH for arbitrary given acceleration and angular 
velocity, van Putten and Eardley use the Lorentz gauge 
condition to obtain second-order wave equations for the 
connection coefficients, van Putten gives a first-order 
form of these equations, involving the Riemann tensor as 
well as the connection coefficients, for his numerical im- 
plementation in ID Gowdy wave vacuum spacetime [l6| . 
The van Putten-Eardley formulation allows for complete 
freedom in the choice of spacetime foliation. 

An example of using orthonormal frames in spatial hy- 
persurfaces is the Ashtekar formulation (with subsequent 
modifications) [13, [H, O |2(1 |2l|]. Ashtekar follows a 
3 + 1 split in the sense that his variables are defined rel- 
ative to a hypersurface. However, his description of the 
geometry of the hypersurface is in terms of orthonormal 
triads instead of the metric. He uses complex variables 
to give a compact formalism, requiring the use of reality 
constraints to recover real spacetime. As subsequently 
modified by Yoneda and Shinkai, Ashtekar's formulation 
becomes a FOSH system of PDE's. Shinkai and Yoneda 
[22I I23I ] have published numerical studies using this for- 
mulation in ID planewave vacuum spacetime. 

II. VARIABLES 

Throughout this paper, lower case Greek letters de- 
note spacetime indices (0-3), and lower case Latin letters 
denote only spatial indices (1-3). The letters in the be- 
ginning of the alphabets, (a, (3, 7, 8, e) and (a, b, c, d, e, 
f), denote tetrad indices. Mid-alphabet letters (A, fi, v) 
and (i, j, k, 1) denote coordinate indices. Repeated in- 
dices are summed in all cases. 

For an orthonormal tetrad of basis vector fields, there 
are twenty- four distinct connection coefficients. These 
coefficients are scalar fields under coordinate transfor- 
mations, and are defined as 

r a /3 7 = e Q ■ V 7 ep — — Tparf, (1) 

with V the covariant derivative operator. The Tnii-y axe 
the same as Ricci rotation coefficients (see Wald 24] ). 
They can be written in terms of commutators of the basis 
vectors: 

1 r 

r a /3 7 = -^{(e , [e a , ep\) -q Sl + 

(e d , [e a , e 7 ]) rjgp - (e s , [ep, e 7 ]) T] Sa }. (2) 

Just as in the 3 + 1 formalisms, it is convenient to 
make a space-time split in the tetrad formulation. Here, 
however, the split is relative to the timelike congruence 



3 



defined by the tetrad rather than the constant-t spacelike 
hypersurface. The connection coefficients can then be re- 
labeled as 3D quantities with spatial triad indices (see 
Wahlquist ^ and ERW). To begin, 



K 



ba 



0?j. 



(3) 



where the symmetric part of K\, a is the rate of strain of 
the congruence, and the antisymmetric part, 



(4) 



is the vorticity of the congruence. If fl a — 0, the con- 
gruence is hypersurface orthogonal, and Kb a is the tradi- 
tional extrinsic curvature of the orthogonal hypersurface. 
Note that the sign of K^a here is the same as in ERW 
and WakL but opposite to that of Misner, Thorne, and 
Wheeler |2^. 

The spatial tetrad connection coefficients can be ex- 
pressed more compactly as a two-index quantity defined 
by 



Nk 



2 E-acd Ic 



(5) 



The diagonal components, N±x, N 2 2, and 7V 33 , describe 
the twists of the spatial triads along the 1, 2, or 3 direc- 
tions, respectively. The combinations N a b + Nba of the 
non-diagonal components represent gravitational wave 
degrees of freedom. It is sometimes convenient to rep- 
resent the antisymmetric part of Nb a by its own symbol 



n a = £ £ abc Nb c . 

The acceleration of the congruence is 



(0) 



(7) 



and the angular velocity of the spacelike triads relative 
to Fermi propagated axes is 



1 



Sabc T 



cbO- 



(8) 



There are nine K^ a and nine iVf, , giving eighteen pri- 
mary variables to be evolved. The three a a and three u a 
are gauge quantities, which, in this paper, are evolved by 
either the Nester or Lorentz gauge, or kept fixed. 

Numerical calculations are performed with a particular 
choice of coordinates x^ . Hyperbolic evolution consists in 
calculating variables on the spacelike hypersurfaces char- 
acterized by a; = t 2 from an initial state specified by 
values of the variables on an earlier spacelike hypersur- 
face x° = ti. Let A£J denote the transformation matrix 
between coordinate basis vectors and the tetrad basis vec- 
tors: 



(9) 



From Eq. © , we obtain directional derivatives along the 
tetrad directions in terms of the partial derivatives along 
the coordinate directions, 



(10) 



The are the coordinate components of the tetrad basis 
vectors. The coordinate metric is constructed from these 
tetrad vector components as = rf 1 ^ \£ X%. 

We find it simpler to use as variables not the sixteen 
A^ , but the following sixteen quantities which completely 
determine the tetrad vector components: nine B k , the co- 
ordinate components of projections of the spatial tetrad 
vectors into the hypersurface; three A a , which measure 
the tilts of the spatial tetrad vectors relative to the hy- 
persurface and are (minus) the tetrad components of the 
3- velocity of the hypersurface frame relative to the tetrad 
frame; the tetrad lapse function a and the three coor- 
dinate components of the tetrad shift vector p k , which 
describe the evolution of the coordinates relative to the 
tetrad congruence. In this paper we evolve the B k and 
the A a as dynamic variables, but take a and the (3 k to 
be fixed functions of the coordinates. Eventually we may 
want to expand the hyperbolic system to include dynamic 
equations for the lapse and the shift. 

In a 3 -I- 1 formalism, the lapse function N is the rate 
of change of proper time with respect to coordinate time 
along the hypersurface normal, and the shift vector N k is 
the rate of displacement of the spatial coordinates with 
respect to coordinate time relative to the hypersurface 
normal, such that the coordinate velocity of the normal 
worldline dx k /dt = —N k . Our tetrad lapse function a is 
the rate of change of proper time with respect to coor- 
dinate time along the tetrad congruence, and our tetrad 
shift vector (3 k is the rate of displacement of the spatial 
coordinates relative to the tetrad congruence worldlincs 
per unit coordinate time. For simple comoving coor- 
dinates, the spatial coordinates are constant along the 
congruence, and (3 k = 0. However, comoving coordi- 
nates are not desirable in black hole calculations, since 
with a finite acceleration tetrad worldlines will be con- 
tinuously advected inward across the event horizon. The 
tetrad lapse must be chosen so that the constant-time hy- 
persurfaces remain spacelike, which is equivalent to the 
condition A a A a < 1. 

The projection of e a into the hypersurface is done 
along the congruence worldlincs, such that 



e a = A a e + B a . 



(11) 



Fig. [T] shows how these various vectors are related from 
the point of view of a frame at rest with respect to the 
hypersurface. 

The directional derivative along a spatial tetrad direc- 
tion is 



D„ 



A a D 



a dx k • 



(12) 



4 



III. THE TETRAD EQUATIONS 
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FIG. 1: Decomposition of e a into a vector tangent to the 
hypersurface, B a , and a vector parallel to the congruence, 
A a eo. Displacements of the spatial coordinates, x k , relative 
to the congruence worldline equal (3 k dt, where (3 k is our tetrad 
shift vector. Displacements of x k relative to the hypersurface 
normal n are N k dt, where N is the 3 + 1 shift vector. 



and the directional derivative along the congruence is 

(13) 



a \dt r dx k ) 
Plugging Eq. (JT^ into Eq. O > we g et 



d_ 

dx k 



Aad_ 

a dt 



(14) 



Comparing Eqs. (|T3|) and (fH|) with Eq. (|TU)l . we can read 
off the tetrad vector components. They are 



A'' 









' l_ 

a 




\ k 




- a 



(B k a 



a 



(15) 



The 3 + 1 lapse function, N, and shift vector compo- 
nents, N k , can be expressed in terms of a, (3 k , A a , and 
B k by constructing g tl and g tk from the tetrad vector 
components using g^ — r\ a ^ Xg. With the relations 



tk 



we get 



N 



a 



VI - A a A a ' 



N 



a A a B k 



1 A a A a 



(17) 



(18) 



The tetrad lapse is smaller than the 3+1 lapse due to the 
time dilation of the tetrad observer in the rest frame of 
the hypersurface. 



In this section, we present the evolution and constraint 
equations for the connection coefficients defined in Sec. 
IITI in terms of directional derivative operators along the 
tetrad directions. The structure of the equations in this 
form is deceptively simple. One has to keep in mind that 
partial differential equations are solved in numerical rela- 
tivity, so the directional derivatives have to be expanded 
into partial derivatives using Eqs. (|13|) and l|14|). In terms 
of directional derivatives, the equations are quasi-FOSH. 
We will present the necessary steps to obtain the true 
hyperbolic form in Sec. IIVI Furthermore, we call the 
"constraint" equations presented in this section quasi- 
constraint equations, as is done in 0, because although 
they contain only spatial directional derivatives, these 
spatial directionals contain time partials (see Eq. I|14|l). 
The details of converting the quasi-constraint equations 
into true constraint equations are given in Appendix lAl 

In both the tetrad and 3 + 1 approaches, there is an 
ordering ambiguity inherent in the derivation of a first or- 
der hyperbolic evolution system from the Einstein equa- 
tions. In first order 3 + 1 formalisms, the spatial deriva- 
tives of the metric are independent variables, Djki (which 
equal ^djhki, where h^i is the spatial metric). Index 
re-ordering of the Riemann tensor is an interchange of 
spatial partial derivatives such that 



diDjki = djD 



ikl ■ 



(19) 



Exploiting this freedom leads to a wide variety of hyper- 
bolic formulations 27]. The standard energy and mo- 
mentum constraint equations can also be added to the 
evolution system to get additional hyperbolic or non- 
hyperbolic structures. In the context of a tetrad ap- 
proach, the connection coefficients are given in terms of 
the commutation relations, Eq. @. To derive integrabil- 
ity conditions from Eq. (0) is messy. We find it is much 
easier to approach the question of ordering ambiguity 
from the symmetries of the Riemann tensor, especially 
since these symmetries are explicit when the indices are 
all up or down, and it is trivial to raise and lower indices 
with the Minkowski metric. We do not consider a whole 
range of schemes as do Kidder, Scheel, and Teukolsky 
[271 , but rather a particular scheme which leads to a sim- 
pler hyperbolic structure than the others. To focus on 
the quasi-FOSH structure of the equations, we present 
only the principal terms here. The second order source 
terms are given in Appendix 1X1 

The Riemann tensor projected onto a tetrad is 

Raf3jS — e a ■ (V 7 V,5 — V«5V 7 — V[ e ^,ej]) e fl 

= D~ t T a fis — D r Q) g 7 + T ae j ^ e ps ~ 

F ae 5 r £ ) g 7 + F Qj 3 e (T^^g — T^), 

(20) 

where D a represents directional derivatives along the 
tetrad directions. The antisymmetry of R a /3j5 on the 
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second pair of indices is explicit in Eq. (|20|l . The anti- 
symmetry on the first pair is also trivial because of the 
antisymmetry of r Q( g 7 on the first two indices. However, 
the Riemann identities 

RafifS = R~fSap, (21) 



RafSjS + RaSPj + RajS0 = (22) 

lead to new Riemann constraints, which we exploit in the 
following. 

First we derive all the possible quasi-constraint equa- 
tions, noting the use of Eq. 14.7 from Misner, Thorne, 
and Wheeler. The energy quasi-constraint equation is 

Goo = #1212 + #2323 + #3131 ~ 2 D a n a . (23) 

An analogous quasi-constraint equation for f2 a is derived 
using the cyclic identity, Eq. (|22() . on Oabc: 

= #0213 + #0321 + #0132 ~ 2 D a Q a . (24) 

The momentum quasi-constraint equations are obtained 
from Goa = 0: 

Gqi = i?0212 + #0313 

~ -D x (K 22 + K 33 ) + D 2 K 12 + D 3 K13, (25) 

G02 = #0121 + #0323 

~ -D 2 (An + A 33 ) + £>i A 21 + D 3 A 23 , (26) 

Go3 = #0131 + #0232 

~-D 3 (K n +K 22 ) + D 1 K 31 + D 2 K 32 . (27) 

Similar quasi-constraint equations involving N ab are de- 
rived solely from Eq. (|21|l applied to spatial indices of the 
Riemann tensor: 

#1213 - #1312 ~ Dx (N 22 +N 33 ) - D 2 Nx2 - D 3 N 13 , (28) 



#2321 - #2123 ~ D 2 (N n +N 33 ) - Dx N21 - D 3 N 23 , (29) 



R 3 i32-R 3 23i~D 3 {N 11 +N 22 )-D 1 N 31 ~D 2 N 32 . (30) 

We now to turn to evolution equations. By using 
Eq. H21(l in the momentum quasi-constraint equations, we 
calculate evolution equations for the non-diagonal com- 
ponents of N ab . For example, interchanging the first and 
second pairs of indices in Eq. (|27|l as shown below gives 
evolution equations for N 2 x and N\ 2 : 

G03 = #0131 + #3202 

~-D N 21 -Diux-DsKu+DtKn, (31) 

Go3 = #3101 + #0232 

- D N 12 + D 1 u 2 + D 2 K 32 - D 3 K 22 . (32) 



Evolution equations for the diagonal components of N a b 
are obtained solely from applications of Eq. (|21|l : 

#2301 - #0123 ~ -Do JVu +Di (j! + D 2 K 31 - D 3 K 21 , (33) 
#3102 - #0231 - D N 22 + D 2 lu 2 + D 3 Ka-Di K 32 , (34) 



#1203 "#0312 - DoN 33 + D 3 u 3 + D x K 23 -D 2 K l3 . (35) 

Evolution equations for K a b are obtained from the Ein- 
stein equations. As is done elsewhere (see for example 
|28j). a multiple of the energy quasi-constraint equation 
times the spatial metric can be subtracted from the K a b 
evolution equations. In the formulation we are present- 
ing, this amounts to subtracting a multiple of Eq. Ij23(l 
from the diagonal spatial components of the Ricci tensor, 
since our metric is just 8\. The same type of procedure 
in 3 + 1 formulations affects the evolution equations for 
the non-diagonal as well as the diagonal components of 
Kab, since the spatial metric there is, in general, non- 
diagonal. The number we choose to multiply the energy 
quasi-constraint equation by is 1, to bring the evolution 
equations for the diagonal K a b into the same form as 
those for the non-diagonal components. With the aid of 
Eq. H21|) . we derive sample evolution equations for K a b 
below: 

#12 = "#1020 + #1323 = "#1020 + #2313 

~DoK 21 -D 2 a 1 -D 3 N n +D 1 N 31) (36) 

#12 = —#2010 + #1323 

~ D K 12 -D ia2 -D 2 N 32 + D 3 N 22 , (37) 

#11 " Goo = "#0101 " #2323 

~ -Do K tl - Dx a x - D 2 N 31 + D 3 N 21 . (38) 

If the gauge quantities a a and Lu a are fixed functions 
of time, then the quasi-FOSH structure is now complete 
and is represented by the following system of eighteen 
equations for the nine K a \, and nine N a b % . 

D K ab ~ D a a b - Eacd D c N db = S_K abl (39) 

D N ab + D a uj b + e acd D c K db = S_N abl (40) 

where S_K ab and S_N ab are source terms quadratic in 
the connection coefficient variables. (See Appendix iBl) 

Alternatively, one can implement a dynamic gauge. 
Both the Nester gauge used by ERW and EW and the 
Lorentz gauge used by van Putten and Eardley result in 
evolution equations for a a and io a which , when added to 
Eqs. ||s2]J and |@DJ|, form a quasi-FOSH system. 

The Nester gauge conditions in Q are defined for an 
arbitrary number of dimensions. In 4D spacetime, the 
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Nester conditions state that two 1-forms, q and q, whose 
tetrad components are 



(41) 



are closed. This implies vanishing exterior derivatives so 
that, in tetrad component form, the Nester conditions 
are 



V Q qp - q a = 0, 



(42) 



V a «/j- Vpqa = 0. (43) 
The Nester conditions result in six evolution equations, 

D a b - D c K cb = -(w c - 2 Sl c ) N cb + (TrN) uj b , 

D uj b + D c N cb = -(oj c - 2 n c ) K cb - {TrN) a b , (44) 

and six quasi-constraint equations, 

e abc D b (o c - 2 n c ) = 2 (TrJf) fi„ - (TrN) (a a - 2 n a )+ 

(a d - 2 ra d ) AT da , 
e abc D b (lu c - 2 O c ) = —(TrN) uj a + (uo d - 2 fi d ) JV*,. 

(45) 

The Lorentz gauge condition in 

7] 

is V 7 w 7a/ g = 0, 

where cj 7C( ^ = r a/ g 7 . For fixed a and /?, the w 7Ct/ 2_are 
the components of connection 1-forms (see Wald (2J|)- 
Expressed in terms of r Q( g 7 , the Lorentz gauge condition 
is 



D5TJ+TJT\ 5 = Q. (46) 

The Lorentz gauge results in six evolution equations: 

Do a b - D c K cb = (a c — 2 n c ) K cb - (Tr-if) a b , 
D uj b + D c N cb = -(a c -2 n c ) N cb - (TrK) w b , (47) 

and no additional constraint equations. Note that the 
principal terms in Eq. (|47|) are identical to those for the 
Nester gauge evolution equations, Eq. 1)44(1. 

The evolution equations expressed in terms of direc- 
tional derivatives, Eqs. 13911 . (|40fl . and (with a dynamic 
gauge) 144|) or (|47ll . can be written in a condensed nota- 
tion, 



Doq + M a D a q=S, 



(48) 



where q is a vector of the twenty-four variables 
N ab , K ab , a a , and u> a . Furthermore, M a are three sparse 
24 x 24 matrices whose only nonzero elements are ±1, and 
S is a vector of source terms. If one orders the variables 
so that 

q = (N 11 ,N 2 i,N 31 ,a 1 ,K 11 ,K 2 i,K 31 ,uj 1 , 

Ni2,N 2 2,N 3 2,a2,Ki2,K 2 2, K 32 ,uj 2 , 

N 13 , iV 23 , 7V 33 , a 3 , if 13 , K 23 ,K 33 ,u> 3 ), (49) 



then the M a matrices have a simple block diagonal struc- 
ture. Each is composed of three identical 8x8 blocks, 
respectively: 
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(50) 



(51) 



(52) 



IV. 



HYPERBOLIC STRUCTURE OF THE 
COORDINATE EQUATIONS 



The beautiful quasi-FOSH tetrad formulation with 
constant coefficients given in Eqs. 1)48(1 to 1(52(1 is not so 
simple, or so beautiful, when expressed in terms of coor- 
dinates. However, we proceed to show that the coordi- 
nate equations are still quite manageable and are, indeed, 
symmetrizable hyperbolic. 

Substituting Eqs. U3J) and (JTIJ into Eq. (gHJ, we ex- 
press the tetrad evolution equations in terms of partial 
derivatives along coordinate directions: 



where 



TD a q + M a B h a ^ = S, 



[I + M a A a ] , 



(53) 



(54) 



and J is the identity matrix. For clarity of notation, we 
let D represent the partials in Eq. 1(13(1 . The T matrix 
has a block diagonal structure composed of three identi- 
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cal 8x8 entries which are 



-*- block 
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-A 3 


A 2 


Ax' 


(J 
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A 3 





-Ax 


A 2 
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-A 2 


Ax 





A 3 
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-Ai 


-A 2 


-A 3 








A 3 


-A 2 




1 











-A 3 





Al 


-A 2 





1 








A 2 


-Ax 





-A 3 








1 





Ax 


A 2 


A, 














1 



(55) 



We now multiply Eq. H53|) by T 1 to give 



D q+C a B k 



d q 

<9x fc 



T _1 S, 



where 



C a ee T~ l M a . 



(56) 



(57) 



T 1 is straightforward to calculate, since it is a block 
diagonal matrix with each block equal to the inverse of 
T b iock in Eq. (J5SJ, where 



rp-L 

block 


1 - 


A A X 
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-A 2 


-A^ 





1 





-A 3 





Ai 


-A 2 
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A 2 


-Ax 





-A 3 
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Ai 


A 2 


A 3 








-A 3 


A 2 Ai 
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A 3 





-Ax A 2 
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-A 2 


Ax 


A 3 
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-A 2 


-A 3 











1 



(58) 



The system given in Eq. (|56() is hyperbolic according to 
the definition given in if any linear combination b a C a 
of C , C 2 , and C 3 can be diagonalized with a complete 
set of eigenvectors and real eigenvalues. The linear com- 
bination for propagation in the x k coordinate direction 
has b a = B k . Solving the eigensystem of the combined 
8x8 matrix in Mathematica gives for the eigenvalues 



-b A ± ^/b- b (b x A) ■ (b x A) 
I- A- A 



(59) 



where b and A are the 3-vectors in the spatial orthonor- 
mal frame with components b a and A a , respectively. The 
dot and cross products are the standard vector opera- 
tions, and the upper/ lower signs on the square root are 
for left/ right- propagating (relative to b) modes. The 
eigenvalues are real as long as A ■ A — A a A a < 1. The 
eigenvectors (given in Appendix [5| form a complete set. 
The lapse and shift hidden in the Do operator modify 
the eigenvalues in a trivial way (multiply by the lapse 
a, then subtract the shift /3 fc ), but have no effect on the 
eigenvectors. 

For Eq. (|56|) to also be symmetrizable by the definition 
of Lindblom and Scheel |12( , a positive definite symmetric 



matrix must be found which multiplies the C a matrices 
to give symmetric matrices. The obvious candidate for 
such a symmetrizer is T, since T C a = M a by Eq. J37J), 
and the M a matrices are symmetric. T is real and sym- 
metric so a necessary and sufficient condition for it to 
be positive definite is that all its eigenvalues are positive 
|29j . This requires that A a A a < 1. 

Note that our saying the system is hyperbolic is con- 
tingent on the evolution of the B k components being hy- 
perbolic. This is discussed further in Sec. 



V. EVOLUTION AND CONSTRAINT 
EQUATIONS FOR A a AND B k a 

The commutator of the orthonormal basis vectors is 

[e a , ep] =\7 a e fj - Vp e a = (T 1 f3a - r 7 Q/3 ) e 7 . (60) 

Expressed in terms of tetrad components and partial 
derivatives, 



Expanding Eq. I|61|l and collecting terms gives 



(61) 



d\ k 




8% 




dt 


dx k 


d\ p 






x* aA « 
~ ~dt 


dt 


dx k 



d\ k 



d 



x 



dx k J dx k 
dX l a \ d_ 
8x k ) dt' 



(62) 



Doing the same in Eq. (|60|l . we get 



[( r V - T °al3) A + ( TC (3a - ^afj) ^c] ~Q^k + 



0a 



a/3) 



d_ 



(63) 



Set Eqs. (|6^|l and l|63|l equal, and let the index a = 
and the index (3 — a. Simplifying, we obtain evolution 
equations for B k and A a : 



k B a 0p n 1(8 \ k 

= -e abc Lu b B k -K ac B k , (64) 



D A a — a a - e a bc uJb A c - K ac A c 



(65) 

a ox 1 



where Cp in Eq. (|64|) is the Lie derivative. For fixed lapse 
and shift, the evolution of both the B k and the A a is just 
advection along the congruence worldlines and trivially 
hyperbolic. 

The congruence can always be chosen to be orthogonal 
to the initial hypersurface, so A a = initially. However, 
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A a will not remain zero during the subsequent evolution 
unless the condition 

Q 

a ~dx ll ° 9a ^ Qa ^ 

is satisfied at all times. For either the fixed or the dy- 
namic gauge conditions on a a considered here, the evolu- 
tion of B l a and/or a a is inconsistent with Eq. I|66|l . except 
possibly for very special initial conditions. 

Repeating the same process as above, but with the 
index a — a and the index = b, we obtain constraint 
equations for B k and A a : 

0B k 

Scab B a -Q-j- = N dc B d - 

(TrN) B k c + e cab A a (e bdf uo d B) + K bd B k d ), (67) 
£cab B l a -^-j = 2 fi c + N dc A d — 

(TrN) A c + e cab A a (-a b + e bdf Lo d A f + K bd A d ). (68) 

Note that if A a = 0, the constraint of Eq. I|68|l is satisfied 
automatically, and Eq. I|67ll can be used to calculate all 
of the N ab from the B k . 

VI. THE INITIAL VALUE PROBLEM 

Initial conditions for the variables must be chosen so 
all relevant constraints are satisfied on the initial hyper- 
surface. This is most easily accomplished if the initial 
tetrad is oriented so the tetrad congruence is orthogonal 
to the initial hypersurface. Then, fl a = and A a = 
initially. The initial B k are the components of an or- 
thonormal triad of vectors lying in the initial hypersur- 
face, related to the inverse of the spatial metric of the 
hypersurface by h kl = B k B l a . One way to construct 
consistent initial conditions is to solve the initial value 
problem using standard 3 + 1 methods for the spatial 
metric and the extrinsic curvature. Construct orthonor- 
mal triad fields by a Gramm-Schmidt orthogonalization 
procedure, orienting the B\ vector along the x 1 coordi- 
nate direction, and Bi in the x 1 — x 2 plane, for instance. 
Find the N ab from the commutators of the B a 3-vectors 
using Eq. (|67|l . The K ab are simply the projections of 
the coordinate components of the extrinsic curvature as 
found from the 3 + 1 initial value problem along the or- 
thonormal triad vectors. This procedure is guaranteed 
to give consistent initial conditions for the tetrad vectors, 
as long as the spatial coordinates are basically Cartesian, 
ie., there are no spatial coordinate singularities. The ini- 
tial acceleration and angular velocity of the tetrad are 
arbitrary, except in the context of the Nester gauge. In 
this case, the initial angular velocity uj a and the initial 
a a — 2 n a must have vanishing exterior derivatives. How- 
ever, there is no guarantee that there are not large twists 
of the initial triad vectors, possibly leading to large gauge 



transients in the context of one of the dynamical gauge 
conditions. 

A more elegant and, likely, a better-behaved choice for 
the initial spatial triad is to require that it satisfy the 3D 
Nester gauge conditions in the hypersurface. These con- 
ditions are that the 3D one- form q ai whose triad compo- 
nents are T a b b , has zero exterior derivative, and that the 
trace of the N ab matrix vanish, ie. Nu + N22 + N33 = 0. 
For simple topologies, the first condition is equivalent to 
the condition that 2 n a = e abc N bc be the gradient of a 
scalar. Finding a solution for the B k which satisfies these 
conditions is, in general, a non-trivial elliptic problem 
(see [3(j and HU). However, the situation is much sim- 
pler for conformally flat 3-geometries. Taking Cartesian 
basis vectors in the conformal geometry, the conformal 
N ab are zero. If the conformal factor which generates 
the physical metric is e 4 ^, the conformal transformation 
simply rescales the B k by a factor e -2 ^. The physical n a 
equal —2 D a ip (the gradient of a scalar), and the sym- 
metric part of N ab is still zero. The 3D Nester condition 
is satisfied. 



VII. DISCUSSION 

Although the emphasis of the tetrad formalisms in the 
literature is on evolving spacetimes with physically de- 
fined flows, we think a tetrad formalism may be useful for 
vacuum numerical calculations of black hole spacetimes. 
The tetrad formulation we have presented, based on that 
of Estabrook, Robinson, and Wahlquist 3], allows con- 
trol of the evolution of the timelike congruence through 
either dynamic or fixed gauge conditions on a a and uj a . 
The lapse and shift we have introduced allow for a com- 
pletely arbitrary evolution of the coordinates with min- 
imal complication of the formalism. Furthermore, since 
the variables evolved are defined relative to the orthonor- 
mal frames, the metric is the Minkowski metric and there 
are no nonlinearities in the equations associated with the 
inverse metric. The system of equations based on coordi- 
nate derivatives is symmetrizable hyperbolic, though ad- 
mittedly more complicated than the quasi-FOSH system 
based on directional derivatives. Finally, the variables 
N ab , K ab , a a and uj a are all scalars, so derivatives of the 
shift only appear in the evolution equations for B k . We 
have successfully implemented this tetrad formalism for 
ID colliding gravitational plane waves [3i|, where it re- 
sults in substantially better accuracy and stability com- 
pared with our calculations in the 3 + 1 framework • 

There are still important issues to be worked out in 
order to apply this tetrad formalism to 3D black hole 
codes. First, one must deal with the complications that 
arise when evolving the congruence as well as the hy- 
persurface. The congruence stays timelike as long as the 
acceleration is bounded; however, in order for the the hy- 
persurfaces to remain spacelike, the condition A a A a < 1 
must be satisfied at all times. If A a A a = 1, the sys- 
tem breaks down completely. The hypersurfaces become 
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null, causing the 3 + 1 lapse to blow up. In addition, the 
coordinate equations become singular (see Appendix IB")l . 
A a A a can be kept small by an adjustment of the lapse 
as the calculation proceeds. Such a resetting of the lapse 
invalidates theorems which bound the growth of the solu- 
tion based on the symmetrizable hyperbolic structure of 
the system, but does not affect other advantages of the 
hyperbolic formulation, such as dealing with boundary 
conditions. 

Second, attaining a stationary solution at late times 
for general 3D black hole numerical calculations may 
be problematic in a tetrad formulation. Such solutions 
are potentially desirable for long-term stability. With 
bounded acceleration, congruence worldlines just outside 
the apparent horizon will necessarily cross the horizon 
and be trapped by the black hole. Worldlines further out 
may actually accelerate outward relative to a static or 
stationary observer. While a shift adjusted just right may 
keep grid stretching under control in spite of this, the ad- 
justment will have to be time-dependent except for very 
special initial conditions on the congruence. Knowing 
what these special initial conditions are requires know- 
ing the whole future solution ahead of time, and thus is 
not practical. 

Both of these issues will probably best be addressed by 
extending the symmetrizable hyperbolic system to allow 
for a dynamic lapse, to keep A a A a small, and a dynamic 
shift, to control grid-stretching. 

APPENDIX A: SOURCE TERMS OF TETRAD 
EVOLUTION AND QUASI-CONSTRAINT 
EQUATIONS 

The tetrad energy quasi-constraint equation, derived 
from Goo = 0, is 

2D a n a = -2 io a n a + N cd N cd + 

\{K cd K dc - (TrK) 2 - N cd N dc - (TrN) 2 ). (Al) 

The tetrad momentum quasi-constraint equations, ob- 
tained from Goa = 0, are 

Di (if 22 + K33) - D 2 K u - D 3 K 13 
= 2 eibc a b Vl c + sue K bd N dc - 2 K lc n c , (A2) 

D 2 (K n + K 33 ) - D 1 K 21 - D 3 K 23 
= 2 e 2 bc a b £! c + e 2 b c K bd N dc - 2 K 2c n c , (A3) 

D 3 (K n + K 22 ) - D 1 K 31 - D 2 K 32 
= 2 e 3bc a b VL C + e 3bc K bd N dc - 2 K 3c n c . (A4) 

Analogous quasi-constraint equations for f2 a and N ab are 
2 D a il a = 2 a a n a + 4 n a Ct a , (A5) 



D l (N 22 + N 33 ) - D 2 N 12 - D 3 N 13 = -2 e lbc u b 

2 K lc n c - e lbc N lb N cl + 2 m {N 22 + N 33 ), (A6) 



D 2 {N n + N 33 ) - D l N 21 - D 3 N 23 = -2 e 2bc oo b H c + 
2 K 2c n c - e 2bc N 2b N c2 + 2 n 2 (N n + N 33 ), (A7) 



D 3 (N n + N 22 ) - D l N31 - D 2 N 32 = -2 e 3bc u b « c + 
2 K 3c fl c - e 3bc N 3b N c3 + 2 n 3 (JV n + N 22 ). (A8) 

To obtain the true constraint equations from these 
quasi-constraint equations, Eq. (|12() must be substituted 
for each appearance of D a . Then Eq. Ij56(l must be used 
to eliminate the terms involving Dq that have been in- 
troduced. The true constraints are only needed to check 
accuracy as the solution evolves, once consistent initial 
data is obtained. 

The sources for K a b and N a b evolution Eqs. (|39|l and 
iflH)) are 

S_K a b = a a ab + Ebcd (—N ac a-d + K ac uj d )+ 

£acd K cb LU d + - Eadf Sbce (K dc Kf e - N dc Nf e ) + 

(TrN) N ab - N ca N cb - K ac K cb + 2 uj b Cl a , (A9) 
S_N ab = -a a uj b + e b cd (K ac a d + N ac u d )+ 

£acd N cb U! d + E ad f £bce N dc Kf e — 

(TrN) K ab + N ca K cb - K ac N cb + 2 a b Q a - (A10) 

APPENDIX B: DETAILS OF COORDINATE 
EQUATIONS 

In this appendix, we expand Eq. I|5()l) . Because of the 
block diagonal structure of the C a and M a matrices, it 
is only necessary to work with one set of eight variables. 
We will let the free index d = 1, 2, or 3 (for the first, 
second, and third set of eight variables in the vector q), 
ie. 

qd = (N ld , N 2d , N 3d , a d , K ld , K 2d , K 3d , u d ). (Bl) 

So as to include a non-zero shift with minimal notation, 
we use Dq to represent the partial derivatives in Eq. JTSj. 
Then we can write Eq. (p)FJ|) as 

D q_d+ 

[C 'liock Bi + C\ lock B\ + C\ lock Sf] — —j: = S d , (B2) 
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where 
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(B3) 
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(B4) 



The expressions for S_Nid, S_N 2 d, S_N 3 d, S_K\d, 
S_K 2d , S_ K 3d in Eq. ljB6)l are obtained from Eqs. (|X9)l 
and (|A10(I . Those for S_a,d and S_uJd are from either 
Eqs. (gU) or Eqs. (|4"7|) . 

The eight eigenvectors of the arbitrary linear combi- 
nation b a C biock consist of four pairs of left and right- 
propagating modes. Each pair only involves one of 
the four variables Kid, K 2 d, K 3 d, ^d- For propaga- 
tion in the x k coordinate direction, b a — B\. As in 
Sec. I1VI we simplify notation by using b and A to 
denote the 3-vectors in the spatial orthonormal frame 
with components b a and A a . Then the u>d eigenvec- 
tors, normalized so wj = b ■ b = \b\ 2 , have N b d — 

[b x (6 x A)] b ± b b \J\b\ 2 -\bx A\ 2 , and a d = 0. 

The K a d eigenvectors, normalized so K a d — |b| 2 , have 

N bd = e abc ([bx(bxA)] c ± b c ^J\b\ 2 -\bx A\ 2 ^j and 

a d = - ([b x (b x A)] a ±b a ^J \b\ 2 -\bx A\ 2 ^j . The 

upper sign on the square root corresponds to the up- 
per sign in Eq. ijSTfl) for the eigenvalues. Note that N a d 
is zero in the eigenvectors. 
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and 



A 3 5_K 2d - 
-A 3 5_ITid 
A 2 S_/r ld - 



A a 

- A 2 S_X 3d 4 
■f Ai S_K 3d 

- Ax S_if 2d 4 



S_N\d - A\ S_ujd ' 
- S_N 2 d - A 2 S_uJd 
S-N 3d - A 3 
S_a d + Ax 5_Xx d + A 2 5_/f 2d + A 3 S_K 3d 
Ax S_a d + S_^xd - A 3 S_N 2d + A 2 S_iV 3o 
A 2 S_a d + S_X 2d + A 3 S_N ld - Ax SJV 3a 
A 3 S_a d + S_K 3d - A 2 S_N ld + A x S_N 2d 
-Ax S_N ld - A 2 S_N 2d - A 3 S_N 3d + S_uj d 



(B5) 



(B6) 
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